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Abstract — A numerical model based on the finite-difference 
time-domain method is developed to simulate fluctuations which 
accompany the dephasing of atomic polarization and the decay of 
excited state's population. This model is based on the Maxwell- 
Bloch equations with c-number stochastic noise terms. We 
successfully apply our method to a numerical simulation of the 
atomic superfluorescence process. This method opens the door to 
further studies of the effects of stochastic noise on light-matter 
interaction and transient processes in complex systems without 
prior knowledge of modes. 

Index Terms — Noise, spontaneous emission, stochastic pro- 
cesses, FDTD methods. Maxwell equations 

I. Introduction 

THE finite-difference time-domain (FDTD) method fT^ has 
been extensively used in solving Maxwell's equations for 
dynamic electromagnetic (EM) fields. The incorporation of 
auxiliary differential equations, such as the rate equations for 
atomic populations |2| and the Bloch equations for the density 
of states of atoms |3|, has lead to comprehensive studies 
of light-matter interaction. Although the FDTD method has 
become a powerful tool in computational electrodynamics, it 
has been applied mostly to classical or semiclassical problems 
without noise. 

Noise plays an important role in light-matter interaction. 
Marcuse solved the rate equations for light intensity and elec- 
tron population including noise terms |T| to illustrate the effect 
of noise on lasing mode dynamics [5J . Gray and Roy extended 
the formulation by adding noise to the field equation in order 
to study the laser line shape f6l. Starting from a microscopic 
Hamiltonian, Kira et al. developed a semiconductor theory 
including spontaneous emission to describe semiconductor 
lasers iT]. While considerable progress has been made, these 
models remain in the modal picture. Knowledge of mode prop- 
erties is required to characterize the noise, making it difficult 
to study complex systems in which the mode information 
is unknown a priori. Without invoking the modal picture, 
Hofmann and Hess obtained the quantum Maxwell-Bloch 
equations including spatiotemporal fluctuations fSl. Although 
it was useful to study spatial and temporal coherence in diode 
lasers, this formalism was based on the assumption that the 
temporal fluctuations of carrier density and photon density 
were statistically independent, which often broke down above 
the lasing threshold. A FDTD simulation of microcavity lasers 
including quantum fluctuations was also done recently f9l. 
This simplified model added white Gaussian noise as a source 
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to the electric field. The noise amplitude depended only on 
the excited state's lifetime. The dephasing process, which was 
much faster than the excited state's population decay, should 
have induced more noise but was neglected. 

Our goal is to develop a FDTD-based numerical method 
to simulate fluctuations in macroscopic systems caused by 
interactions of atoms and photons with reservoirs (heatbaths). 
Such interactions induce temporal decay of photon number, 
atomic polarization and excited state's population, which can 
be described phenomenologically by decay constants. The 
fluctuation-dissipation theorem demands temporal fluctuations 
or noise to accompany these decays. We intend to incorporate 
such noise in a way compatible with the FDTD method, that 
allows one to study the light-matter interaction in complex 
systems without prior knowledge of modes. In a previous work 
ifTO I. we included noise caused by the interaction of light field 
with external reservoir in an open system. In this paper, we 
develop a numerical model to simulate noise caused by the 
interaction of atoms with reservoirs such as lattice vibrations 
and atomic collisions. As an example, we apply the method to 
a numerical simulation of superfluorescence in a macroscopic 
system where the dominant noise is from the atoms rather than 
the light field. 

We start with the Bloch equations for two-level atoms in 
one dimension (ID) where the direction of light propagation 
is along the x-axis. 
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where 17/? = jEz/h is the Rabi frequency, ujq the atomic tran- 
sition frequency, the electric field which is parallel to the 
z-axis, 7 the dipole coupling term. Phenomenological decay 
times due to decoherence T2 and the excited state's lifetime 
Ti (which includes spontaneous emission and non-radiative 
recombination) are appended. In the absence of strong light 
confinement, which holds for macroscopic systems, Ti and T2 
can be considered independent of the local density of states 
(LDOS). Hence, they do not have a dependence on spatial 
location nor frequency. We also include incoherent pumping 
of atoms from level 1 to level 2. The rate is proportional to 
the population in level 1, and can be written as PrPii- Pa'*^ 
represents the steady-state value of when E^ = 0. 

The relations between the Bloch vector and the density 
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matrix are 



Pi = Pl2 + P2I 
P2 =i{Pl2 ~ P21) 
P3 = P22 - Pll- 



(2) 



The total polarization of N atoms in a volume V is ^ 
— (A^/y)|7|pi and inserted into the Maxwell's equations 
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The atom-reservoir interactions not only cause decay of 
the Bloch vector, but also introduce noise according to the 
fluctuation-dissipation theorem. In Section [III we describe the 
model developed to include noise in the Maxwell-Bloch equa- 
tions. The FDTD implementation of this model is presented 
in Section ini] In Section [TVl we simulate atomic superfluores- 
cence and compare the results to previous experimental data 
and quantum-mechanical calculations. 

II. Noise Model 

Starting from the quantum Langevin equation within the 
Markovian approximation, Drummond and Raymer derived 
a set of stochastic c-number differential equations describing 
light propagation and atom-light interaction in the many-atom 
hmit IfTTI . The noise sources in these equations are from both 
the damping and the nonlinearity in the Hamiltonian. The latter 
represents the nonclassical component of noise, giving rise to 
nonclassical statistical behavior Since our primary interests 
lie with classical behavior of macroscopic systems, such as 
superfluorescence and lasing, we neglect the nonclassical noise 
in this paper The amplitude of classical noise accompanying 
the field decay is proportional to -^n, where n is the thermal 
photon number. At room temperature the number of thermal 
photons at visible frequencies {Tllo ^ \ eV) is on the order 
of 10^^^. This can be interpreted in a quantum mechanical 
picture as that most of the time there are no thermal photons 
at visible frequencies in the system. Thus, the noise related to 
field decay is neglected in this paper At higher temperatures 
or longer wavelengths, this noise becomes significant and it 
can be incorporated into the FDTD algorithm following the 
approach we developed in our previous work j 10|. 

The classical noise related to the pumping and decay of the 
atomic density matrix can be expressed as 

ri2 = (6 +«C2)V7pP22 
T2I = (6 - i£,2)y/jpP22 



T22 = i-i\J P22lT\ + PrPw). 
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These noise terms are associated with pi2, p2\, and P22 
respectively. 7p = l/r2 — XjlTi. The terms are real, 
Gaussian, random variables with zero mean and the following 
correlation relation 



where j, fc = 1,2,3. The noise terms and represent 
fluctuations corresponding to decoherence by dephasing, while 
r22 is the fluctuation corresponding to relaxation of and 
pumping to the excited state's population. Only the linear term 
for pump noise is included here, a common first order approxi- 
mation |, 12,| . Furthermore, because we assume Ti <C Ti, pump 
fluctuations are neglected in ri2 and r2i since they are orders 
of magnitude smaller than noise due to dephasing. According 
to Q, the noise terms for the Bloch vector are reduced to real 
variables as 



Ti = 2^i^7pP22 

^2 = -26V>^22 

Ta = 2(3^p22/Ti+PrPll). 



(6) 



They can be added directly to ([T]i. 

In a ID system, the total number of atoms N are split 
equally among AI grid cells, giving the number of atoms per 
cell Ns — N/M. All quantities are defined at each individual 
grid cell, e.g. the term P3{x) is the number of inverted atoms 
in one cell at position x. The number of atoms in each 
cell is assumed to be constant assuring pn + P22 = 0. We 
forcibly keep Ng constant via the relation pn — Ng — P22 and 
only calculate the excited state's population P22{t). The final 
stochastic equations to be solved are 
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In the above equation, the steady-state value of p^ in ^ is 
substituted by p^'^ = Ns{TiPr-l) / [TiPr + l), an expression 
obtained by setting the time derivatives in ([TJ to zero, pn in 
the expression of r22 in (IDi can be replaced by Ng ~ P22- 

III. Numerical Implementation 

The most commonly used method of solving the Maxwell- 
Bloch equations is the "strongly coupled method." With At 
being the time step, E and p are both computed at nAt, 
{n + l)At, etc., while H is computed at (n — 1/2) At, 
(n+ l/2)At, etc. This produces equations with coupled terms 
such as that must be solved by a predictor-corrector 

scheme (as used in |3|) or a fixed-point procedure, both 
of which are computationally inefficient. Therefore, we use 
a weakly coupled method that is easily implemented and 
efficient for ID systems. 

The weakly coupled method was put forth by Bidegaray 
ifTSl . The electric field E^ is computed at times nAt, {n + 
l)At, but p is calculated at (n - l/2)At, [n + l/2)At, 
thereby decoupling those discretized equations and creating 
a simple leap-frog type propagation system for ID. The 
noise terms in (|7]) are present throughout the entirety of 
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the simulation and thus, should be incorporated efficiently. 
After discretization, the terms are correlated according 
to {^j{xu,tm)^kixv,tn)) = {1/ At)SjkSuvSmn, and Can be 
generated quickly with the Marsaglia and Bray modification of 
the Box-Miiller Transformation |14|. Because the noise terms 
contain ^/P22, as seen in (|4]i and (|6]l, we are not able to use the 
weakly coupled scheme to solve for pi , p2 and P22 as precisely 
as possible. Instea d, the a pproximation of using the previous 

time step value \J~(^^^^ is employed. It is valid as long as 
the atomic population is varying slowly. For the simulation 
of superfluorescence in Sec. |IV] the maximum change of P22 
over one time step At is only 0.0007%. 
The discretized equations with noise are 
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where we have defined A = |7|/T4eT'2 and B = |7|cjo/14e- 
These equations are solved to obtain the final FDTD equations 
for E^, Hy, pi, p2 and p22- 

IV. Results and Discussion 

We apply the Maxwell-Bloch equations with noise to a 
FDTD simulation of superfluorescence (SF) and compare 
the results to previous data obtained experimentally lITSl 
and theoretically fT6\. SF is the cooperative radiation of an 
initially inverted but incoherent two-level medium resulting 
from spontaneous buildup of a macroscopic coherent dipole. 
This is an interesting and suitable case to study with our 
method because both spatial propagation of light and noise are 
important. Noise caused by collisional dephasing can seriously 
disturb SF and change the emission character to amplified 
spontaneous emission (ASE). We simulate the transition from 
SF to ASE with increasing dephasing rate, corresponding 
to the experiment by Malcuit et al. on super-oxide ions in 
potassium chloride (KC1:02^) |15|. 

Experimentally the ions inside a cylinder of diameter d 
= 80 pm and length L = 1 mm were excited by a short 
pulse. The total number of excited ions is = 3 x 10^. 
The emission wavelength is A = 629 nm. The Fresnel number 



for the excitation cylinder is F = A/ XL ~ 1, where A is the 
area of the cylinder cross-section. Ti = 76 ns, and T2 was 
varied via temperature change. The "cooperative lifetime" or 
the duration of SF pulse = SttATi/SA^A^ is 2.7 ps. The 
estimated delay time for the SF peak after the excitation pulse 

1 ^ 

Td = Tr ^H2^N) (9) 



is 94 ps. 

Since F ^ 1, the EM modes propagating non-parallel to the 
cylinder axis are not supported 1 17|. Those modes propagating 
along the cylinder axis do not have a strong radial dependence, 
nor are there significant diffraction losses. Thus the system 
can be considered as ID in our FDTD simulation. The grid 
resolution is Aa; = 70 nm and the total running time is 
Tsim — 3 ns. The Courant number S is set to 0.999999. The 
magic time step, 5=1, was seen to cause an instability 
in some cases. The value S — 1/2, however, does not 
propagate the large sudden impulses of the noise accurately. 
Setting S = 1 — 10~^ preserves the accuracy to an acceptable 
degree while eliminating the instability at S* = 1. There is 
some numerical dispersion and reflection from the absorbing 
boundary layer, but the error is of the order 10^^. Ignoring 
non-radiative recombination, the atomic dipole coupling term 
I7I = V3A3neo/87r2ri = 1.1 x IQ-^o C-m. 

The simulation is started with the initial condition of all 
the atoms being excited (p22 = Ns). However, because the 
atomic population and polarization operators do not commute, 
the uncertainty principle demands a nonvanishing variance in 
the initial values of the Bloch vector [17|. This results in a 
tipping angle 9 of the initial Bloch vector away from the top 
of Bloch sphere {pi — 0, p2 — 0, ps = Ng). The value of 9 is 
given by a Gaussian random variable centered at zero with a 
standard deviation 9t = 2/\fWs. Since there is no incoherent 
pumping at i > 0, Pr is set to 0. 

Figure [T] shows the output EM energy at a spatial grid 
point outside the system for four different values of the 
dephasing time T2- When T2 — 100 ps > r^, the cooperative 
emission characteristic of SF is clearly seen in Fig. [TJa). The 
number of atoms that emit cooperatively is estimated to be 
Nc = SttcTiA/SA^L = 3.5x 10* and is known as the Arecchi- 
Courtens cooperation number. Since A^c < ^ = 3 x 10^, 
the SF oscillates in time, with the maximal emission intensity 
at t ~ 170 ps. This behavior agrees well with the previous 
result in ||T6| . For T2 = 33.3 ps < r^, there is enough 
dephasing to disturb the cooperative emission. The emitted 
pulse broadens and the time delay increases, as shown in Fig. 
[TJb). For T2 = 25 ps, a further damping of superfluorescence is 
seen in Fig. [Tfc). As T2 decreases more, the pulse continues 
to broaden but the time delay begins to decrease. When T2 
reaches the critical value ^/rPFd = 15.9 ps, the amount of 
dephasing is sufficient to prevent the occurrence of cooperative 
emission. No macroscopic dipole moment can build up and 
the atoms simply respond to the instantaneous value of the 
radiation field. Hence, SF is replaced by ASE. Figure [TJd) 
plots the ASE pulse for T2 = 14.3 ps. The time delay is 
almost immeasurably small and the emission intensity is very 
noisy. Figure |2] compares the delay times taken from our 
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Fig. 1. Numerical results of the output EM energy from initially-inverted two-level atoms, obtained by FDTD solution of the Maxwell-Bloch equations with 
noise. The left three columns show the output energy for three random realizations. The last column on the right shows the output energy averaged over 30 
random realizations. All insets in the last column magnify the temporal range < t < 1 ns. Dephasing time T2 = 100 ps (first row), 33.3 ps (second row), 
25.0 ps (third row), and 14.3 ps (fourth row). 



FDTD simulations to previous results obtained experimentally 
|[T5l and by full quantum-mechanical theory of SF |16l. 
The excellent agreement validates our FDTD-based numerical 
method. We emphasize that inclusion of the noise terms in ([TJ 
is essential to obtain the correct variation of with T2. As 
found in |15|, the previous approach of modeling the initial 
fluctuations as random tipping angles of the Bloch vector and 
ignoring the noise at later times brings about good agreement 
with experiment only when T2 is large making the amplitude 
of the noise terms in (|7]l small. As the dephasing rate increases, 
fluctuations can no longer be modeled simply as an initial 
noise. 

We have also studied the decoherence process. The am- 
plitude of the Bloch vector ps = \/ pf + P2 ~^ P3 — 
^ iV| + ^pvipn — 4/322P11. In the absence of decoherence, 
P12P21 — P22P11, and pb — Ng. The presence of decoherence 
decreases the off-diagonal terms of the density matrix, thus 
P12P21 < P22P11 and Pb < Ng |18|. We estimate the degree 
of decoherence through the ratio ps/pB, which is plotted in 
Fig. |3]for four different values of T2. Each curve is obtained by 
spatial average of pa and pB over the entire excitation region 
and then ensemble-average over 30 realizations. 

When the dephasing time is large {T2 > Td), a macroscopic 
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Fig. 2. (Color online) Comparison of delay times of emission pulse obtained 
by our numerical simulation (black solid circles) with previous experimental 
data (red crosses) and quantum-mechanical calculation results (blue open 
diamonds) taken from ref. |16|. Numerical delay times are obtained from 
the emission pulses averaged over 30 realizations. 
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Fig. 3. (Color online) The ratio ps/ Pb as a function of time for T2 - 100 
ps (red dotted line), 33.3 ps (blue dashed line), 25.0 ps (green solid line), and 
14.3 ps (black dash-dotted line). 

dipole moment is spontaneously formed. The enhanced radia- 
tive decay rate results in quick depletion of the population 
inversion p^. Despite T2 <C Ii, the decay of pi and p2 by 
dephasing is overshadowed by the decay of p^ by SF, leading 
to a rapid drop of p^/ ps in time. This behavior is shown 
by the red dotted line in Fig. [3] The non-monotonic decay 
is caused by SF oscillations as can be seen in Fig. [Ua)- The 
oscillatory SF is a result of the number of atoms being greater 
than the Arecchi-Courtens cooperation number {N > Nc). 
The intensity oscillation leads to an oscillation of population 
inversion which is 90 degree out of phase. The local maximum 
of p3 at t = 320 ps (red dotted curve in |3]l occurs just before 
the second peak of intensity at t = 37Qps [Fig. \lla)]. As 
T2 is reduced, the increased amount of decoherence frustrates 
the buildup of a macroscopic dipole moment and reduces the 
radiative decay rate. Consequently, the depletion of population 
inversion is slowed down. It leads to a slower decay of ps/ps 
and the disappearance of damped oscillations. Finally when 
the dephasing time is small enough (T2 < ^jTrTd), the system 
stays in a decoherent state, and pzj Pb remains close to one 
for a very long time. 

V. Conclusion 

We have developed a FDTD algorithm to incorporate 
stochastic noise in macroscopic systems into the Maxwell- 
Bloch equations. Such noise, resulting from atom-reservoir 
interactions, accompanies the dephasing of atomic polarization 
and decay of and pumping to the excited state population. We 
applied our algorithm to a numerical simulation of superfluo- 
rescence in a ID system. The results are in good agreement 
with previous experimental and theoretical studies. Although 
our simulations only include classical noise, nonclassical noise 
may be incorporated as well. Since they consist of nonlinear 
terms 1 1 1], the incorporation of nonclassical fluctuations to the 
FDTD algorithm may be numerically challenging. Given the 
rapid progress in development of various numerical methods 
of including nonlinearity in the Maxwell-Bloch equations |fT9ll . 
II20I . we are optimistic that the quantum noise terms may be 



successfully integrated into our method. Therefore, our FDTD- 
based model can be used for numerical studies of light-matter 
interaction and transient processes in complex systems without 
prior knowledge of modes. 
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